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Using large-scale molecular dynamics simulations of a two-component Lennard- Jones model in three 
dimensions, we show that the late-time dynamics of spinodal decomposition in concentrated binary 
fluids reaches a viscous scaling regime with a growth exponent n — 1, in agreement with experiments 
and a theoretical analysis for viscous growth. 
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The dynamics of phase separation in multicomponcnt 
fluids involves very rich and general phenomena, and has 
therefore been the subject of intensive studies in recent 
years. The dynamics of first-order phase transitions in 
general, besides being of technological importance, is par- 
ticularly interesting because of the emergence of one char- 
acteristic length scale, R(t), during the late times of the 
dynamics. R(t) is related to the average domain size 
of the ordering phase, and displays a simple power-law 
dependence with time, R(t) ~ t n , where n is the growth 
exponent. The presence of one characteristic length scale 
during the late times leads to an interesting dynamical 
scaling behavior, as can be detected from the density- 
fluctuation pair-correlation function G(r,t) = g[r/R(t)], 
or the structure factor, 5(q, t) = R(t) d F(x), where d 
is the spatial dimension and x — Q-^W is the scaled 
wavevector [ jj . 

Whereas the dynamics of phase separation in alloys, 
with conserved order parameter, is quite well understood 
in terms of the Lifshitz-Slyozov theory [Q , and is charac- 
terized by a growth exponent, n = 1/3, independent of 
spatial dimension, volume fraction p| and even the num- 
ber of coexisting phases ||, the dynamics in fluids is a 
more complicated phenomenon due to the coupling of the 
additional velocity field (which is absent in alloys) to the 
ordering field. Consequently various competing effects 
may appear in phase-separating fluids leading to various 
growth exponents depending on the strength of the cou- 
pling between the velocity field and the ordering field, on 
the volume fraction 0-|(|, on the spatial dimension, and 
even on the number of components . 

There is no satisfactory theory for the phase separa- 
tion dynamics in fluids. Thus our understanding of the 
phenomenon is achieved essentially through numerical 
studies and dimensional analysis of the relevant dynam- 
ical model. Using heuristic arguments, Siggia [Q was 
the first to propose that the growth exponent is n = 1 
in phase-separating binary fluids with relatively com- 
parable volume fractions of the two components. This 



growth regime is due to an instability of the tubular do- 
main structure in binary fluids, leading to the transport 
of material from the necks to the bulges. The numeri- 
cal studies of the phenomenon are mainly carried out by 
means of three different methods: Numerical integration 
of the corresponding kinetic phase-field model known as 
model H ||; lattice-Boltzmann (LB) simulations and 
molecular dynamics (MD) simulations Jl0| ], In contrast 
to the first two methods, in a molecular dynamics sim- 
ulation, the hydrodynamic modes arize naturally from 
the microscopic interactions between the molecules sub- 
sequent to a quench into the fluid phase. There has been 
some concerns with regard to the validity of molecular 
dynamics in studying the late-time dynamics of phase 
separation due to the very small time scale involved. It 
should be noted that phase separation in simple fluids is 
naturally a very fast process. Therefore, in order to probe 
the dynamics, experimentalists must perform very shal- 
low quenches, using the advantage of the increased time 
scale due to critical slowing down. In contrast, quenches 
are very deep in a typical molecular dynamics simulation. 

The numerical integration of model H leads to an 
asymptotic growth exponent, n — 1, in agreement with 
Siggia's prediction. LB simulations also find the same 
result M. However, a recent MD simulation on the 
two-component Lennard- Jones potential by Ma et al. pp| ] 
suggests a growth regime with an exponent very close to 
2/3. As we will see later, such an exponent is due to 
inertial effects, and can be calculated from dimensional 
analysis. A more recent model H simulation by Lookman 
et al. |llj] finds that by decreasing the shear viscosity 
of the fluid, a growth exponent of n = 2/3 can be ob- 
served. We are therefore faced with the problem that 
whereas numerical simulation calculations in the case of 
phase separation in alloys agree with the theoretical pre- 
dictions, numerical simulations which are expected to 
most faithfully describe the true dynamics, i.e., molec- 
ular dynamics simulations, are not in agreement with 
theoretical predictions in the case of phase separation 
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in binary fluids. In order to elucidate this apparent dis- 
crepancy between the previous numerical studies and the 
MD simulations of Ma et al, we have carried out a large- 
scale and systematic molecular dynamics simulation of 
the two-component Lennard-Jones model and found re- 
sults which disagree with the MD simulation of Ma et al. 
but are fully consistent with experiments and previous 
model H and lattice-Boltzmann simulations. It is worth 
noting that the present study is the first large-scale MD 
simulation on three-dimensional binary fluids in which 
the viscous regime is observed. 

In our simulation model, we consider N monoatomic 
molecules interacting through the following two- 
component Lennard-Jones potential: 



with oti = 1 if i is an A-molecule, and o.j = 2 if % is a 
B-molecule. In Eq. 
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is the distance separating the 
molecule from the j th molecule, and r c aiOL . is a cutoff 

distance which is equal to 2.5a for 014 = ctj and 2 1 / 6 cr for 
en ^ ctj. 9{x) is the standard Heaviside function. The 
phase diagram of this model, which has recently been cal- 
culated by means of mean field theory and Monte Carlo 
simulation, has a consolute point at T c ~ (4.7±0.2)e for a 
fluid density of p = 0.8a' 3 p2[ . We have performed crit- 
ical quenches at temperatures keT = 2, 3, 3.5, 3.75 and 
4e as well as off-critical quenches at keT = 2e. Notice 
that we have not made quenches to very low tempera- 
tures in order to avoid the solid-gas coexistence region. 
The temperature is controlled by a Nose-Hoover thermo- 
stat [ [l3f , and the Hamilton equations are integrated using 
the Leap-Frog algorithm with a time st ep of At = 0.005r 
where the time scale is r = y pa 2 / 1 e with p being the 
molecular mass. In all of our simulations, the total num- 
ber of molecules is N — 343 000, an order of magnitude 
lar ger than the largest system size considered by Ma et 
al. flO| . Our simulations were performed on an IBM SP2 
parallel machine using 12 processors. Furthermore, a sta- 
tistical average is performed for each quench; 16 runs for 
ksT = 2e and 4 runs for all other quenches. 

We have calculated both the correlation function 
G(r,t) = (<t>(r,t)ct>(0,t)) , where 0(r,i) = \p A (r,t) - 
PB{r,t)]/p is the order parameter and pa and ps are 
the local densities of the two components. We have also 
calculated the structure factor 5(q, t) — (|</>(q, t)\ 2 )/V, 
where 0(q) is the Fourier-mode of the order parame- 
ter and V is the system volume. Both the structure 
factor and the correlation function are then spherically 
averaged. The average domain size is then defined as 
the first zero of the correlation function, Rc(t), and 
as the n th moment of the structure factor, R n (t) — 
2irlfd q S( q ,t)/fd q q n S(q 7 t)} 1/n : 

The time evolution of the pair-correlation function is 
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FIG. 1. Scaled pair-correlation function, g(x), versus the 
scaled distance, x = r/R(t), for a quench at keT = 2e. The 
data shown range from t = 80 to 220Y. The inset shows the 
time evolution of the correlation function from t — 20 to 220r 
in steps of 20r. 



shown in the inset of Fig. 1 for a quench at kBT = 2e. 
The presence of the decaying oscillations in G(r, i) indi- 
cates the occurrence of phase-separated domains which 
are correlated within short distances, due to the conserva- 
tion of the ordering field. The first zero of the correlation 
function increases with time implying a coarsening of the 
domain structure. We have verified that the system has 
reached a dynamical scaling regime by observing the scal- 
ing of the correlation function, shown in Fig. 1 , for times 
larger than about t — 80r. Very good scaling is also 
observed in the structure factor (not shown). The pres- 
ence of a unique length scale in the system at late times 
implies that the width of the interfaces become vanish- 
ingly small compared to the domain size. As a result, 
the structure factor should scale as for large q, 

which is known as Porod's law and is usually observed in 
phase-separating systems at late times. Indeed, we found 
that the structure factor is consistent with Porod's law, 
implying that the phase separation process in our simu- 
lations is well within a dynamical scaling regime. 

Now that we are confident that the systems, we are 
dealing with in our simulations, are safely within a scaling 
regime, we turn to the discussion of the nature of the 
growth law. In Fig. 2, the time dependence of the average 
domain size, as calculated from the various definitions, is 
shown. Notice the linear dependence of R(t) at late times 
indicating that the growth regime should be viscous, in 
agreement with Siggia's prediction However, when 
plotting the data in a double-logarithmic plot, we find 
that the late-time growth exponent is more consistent 
with 2/3, possibly indicating that the observed growth 
regime is inertial, as suggested by the MD simulation of 
Ma et al. [|o|. It should be pointed out, however, that 
the growth law, R(t) = R(0) + at, investigated over a 
finite time range may show a growth exponent which is 
smaller than 1 due to a non-negligible value of R(0) and 
possible other non-algebraic dependences. 

In order to determine the true asymptotic growth law, 
we have to consider the relevant dynamical model and 
analyze our results in light of its implications. The dy- 
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FIG. 2. The average domain size as a function of time 
for a quench at keT = 2e. Ra(t) is the first zero of the 
pair-correlation function, and R\ (t) is calculated from the first 
moment of the structure factor. The two dots indicate the 
typical size of the error bars in the numerical data, and the 
two dotted lines are straight lines. 



namics of phase separation in fluids can be described by 
the so-called model H corresponding to a general- 
ized Cahn-Hilliard equation coupled to the Navier-Stokes 
equation. The appropriate dynamical equations can then 
be written as follows: 



d t (j){r, t)+v V0(r, t) = MV : 



p [d t v(r, t) + (v(r, t) • V) v(r, t)] = V \7 2 v(r, t) 



-Vp(r, t)-0(r,t)V 



(2) 



(3) 



where </>(r, t), v(r, i) and p(r,t) are the local order pa- 
rameter, the velocity field and the pressure field, respec- 
tively. The constants M , p and rj correspond to the order 
parameter mobility, the fluid density and the shear vis- 
cosity, respectively. T is the usual </) 4 free energy func- 
tional 0. The difference between Eq. (||) and the usual 
Cahn-Hilliard equation is the presence of the second term 
on the left-hand side which accounts for the transport 
of the order parameter by the velocity field. Eq. (|3[) is 
different from the usual Navier-Stokes equation by the 
presence of the additional force acting on the fluid due 
to gradients in the chemical potential. 

The set of equations, (2) and (3), is very difficult to 
solve, but one can obtain various growth regimes by 
means of simple dimensional analysis. Here we will limit 
ourselves to three dimensions. At relatively early times, 
but late enough so that the domains are well defined 
and much larger than the interfacial width, the veloc- 
ity field is decoupled from the order parameter leading 
to the usual Lifshitz-Slyozov growth law usually observed 
in alloys, R(t) ~ (Mjt) 1 ^ 3 , where 7 is the interfacial ten- 
sion [|2j. This regime will be referred to as the diffusive 
regime. At later times, the coupling between <p and v can- 
not be neglected, but the inertial term in Eq. (^), can be 
neglected, so that v becomes slaved by </>. One thus ob- 
tains the following growth law, R(t) ~ ("ft/rj), which will 




FIG. 3. (a) The average domain size as a function of 
t v — {'y/rfjt. Data lines from bottom to top correspond to 
kaT = 2, 3, 3.5, 3.75 and 4e respectively. For the sake of 
clarity, data data has been shifted vertically upwards, (b) 
The average domain size as a function of ti = (7/ p) 1 '' 3 ^; 2 '' 3 . 
Data from top to bottom correspond to keT = 2, 3, 3.5, 3.75 
and 4e respectively. 



be associated with a viscous regime, and has been pre- 
dicted by Siggia [Q as a consequence of a necking-down 
instability of the tubular (interconnected) domain struc- 
ture due to the transport of material from the necks to 
the bulges. This regime has been observed in several sim 
pie binary fluids and binary homopolymer blends Jl5 16 
as well as in numerical simulations At even later 

times, the inertial term in the Navier-Stokes equation can 
no longer be neglected, and one finds the growth law of 
the inertial regime, R(t) ~ (7/p) 1/3 i 2/3 §. The two last 
regimes can be observed only for interconnected domain 
structures. For dilute binary solutions, the domains are 
droplet-like, and the domain growth is essentially due to 
their coalescence leading to a growth law, R(t) ~ i 1 / 3 , 
but with a prefactor which is larger than that in the dif- 
fusive regime. The inertial regime has not been observed 
experimentally, but it has been observed in several nu- 
merical simulations in two dimensions |^,|l^,|l8|. Intro- 
ducing the following two time scales, t v = ("f/r])t and 
ti = (7/p) 1 / 3 f 2 / 3 , for the viscous regime and the inertial 
regime respectively, the ^-dependence (i^-dependence) of 
R(t) during the viscous (inertial) regime must be linear 
and independent of the quench depth, except maybe for 
interference with R(0). 

We have therefore calculated, by molecular dynamics 
simulations, the interfacial tension, 7, and the shear vis- 
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cosity, 77, for the various quench temperatures consid- 
ered in the present study. We obtain a shear viscos- 
ity which is practically independent of temperature, and 
equal to 77 — 1.65. However, the interfacial tension is 
found to decrease with temperature, almost linearly, from 
(1.85 ± 0.07)e/cr 2 for k B T = 2e to (0.39 ± 0.13)e/er 2 for 
keT = 4e, since the present model belongs to the Ising 
universality class in d — 3. In Fig. 3(a), R(t) is plotted 
versus t v = {'yt/rf), and in Fig. 3(b), R{t) is plotted versus 
ti = (7/p) 1 / 3 ?; 2 / 3 for all quench temperatures. Although 
the data is almost linear with ti for all temperatures, the 
slope of R(t) versus ti depends strongly on T, whereas 
the slope of R(t) versus t v is independent of temperature. 
This, therefore, strongly indicates that the growth regime 
found in this system cannot be inertial, in contrast to the 
prediction by Ma et al [ fl0| |, but in agreement with the 
other numerical studies and experiments. In our simula- 
tions, dynamical scaling is observed starting from t = 80r 
at the lowest quench temperatures. At higher tempera- 
tures, the scaling regime is delayed to later times. This 
is to be contrasted to the study of Ma et al, in which 
it was found that the scaling regime starts as early as 
20t Of course, the fact that we did not observe 

an inertial regime does not disprove the presence of this 
regime at even later times, as predicted by the scaling 
analysis. The inertial regime has been observed in previ- 
ous numerical studies in two dimensions [ p~T|JTl| ] , and in 
a recent model H simulation in three dimensions fll|| . In 
order to detect such a regime, we must simulate much 
larger systems. 

Another reason, making us even more confident that 
the dynamical regime found in the present study is vis- 
cous, is the value of the prefactor of the growth law 
in terms of t v . Siggia predicted that this prefactor is 
0.6, whereas San Miguel, Grant and Gunton Q find 
that it should be 0.25 from a linear stability analysis of 
the tubular structure. However, a detailed experimen- 
tal study by Guenoun et al. find that the prefactor is 
0.138 ± 0.006 Jl(J. In our simulation, we find the pref- 
actor to be 0.11 ± 0.01, which is very close to the ex- 
perimental value of Guenoun et al., but disagrees with 
the two theoretical predictions which, however, are quite 
crude in nature. The difference between the value of our 
prefactor and that of Guenoun et al. might be due to 
the finite size of our systems, leading to a cutoff of the 
long-range hydrodynamic modes. Indeed, one expects 
that this prefactor decreases linearly with 1/L from its 
thermodynamics-limit value [fl!^l . Moreover, we found 
that the prefactor of t v is independent of volume fraction 
for quenches at volume fractions around 0.5. However, 
for volume fractions smaller than about 0.3, we found 
a growth exponent consistent with 1/3. We should no- 
tice that recently, Nikolayev et al. have predicted that a 
sharp transition from the viscous growth to coalescence- 
dominated growth occurs at a volume fraction around 
0.3 f§. 

In conclusion, we have performed a large-scale sys- 
tematic molecular dynamics study of the phase separa- 



tion dynamics in binary fluids in three dimensions which 
faithfully accounts for hydrodynamic modes. During late 
times, the system reaches a dynamical scaling regime dur- 
ing which the average domain size grows linearly with 
time in agreement with Siggia's prediction, previous nu- 
merical integration of model H, and lattice-Boltzmann 
simulations. The discrepancy with a previous molecular 
dynamics study has been clarified. 
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